function logistic_error

%  error obtained using various DE solvers for the problem
%  y' = ry(1-y)  with   y(0) = y0      '

% clear all previous variables and plots
clear *
clf

% set parameters for problem
r=10;
t0=0; y0=0.1;
tmax=1;

% loop that increases M 
for im=1:4

	M=10^im;
	points(im)=M;
	t=linspace(t0,tmax,M+1);
	h=t(2)-t(1);

	% exact solution
	y_exact=exact(t,y0,M+1);

	% euler
	y_euler=euler('de_f',t,y0,h,M+1);
	%e_euler(im)=norm(y_exact-y_euler,inf);
	e_euler(im)=abs(y_exact(M+1)-y_euler(M+1));

	% backward euler
	y_beuler=beuler('de_f',t,y0,h,M+1);
	%e_beuler(im)=norm(y_exact-y_beuler,inf);
	e_beuler(im)=abs(y_exact(M+1)-y_beuler(M+1));
   	
	%  trap
	y_trap=trap('de_f',t,y0,h,M+1);
	%e_trap(im)=norm(y_exact-y_trap,inf);
	e_trap(im)=abs(y_exact(M+1)-y_trap(M+1));

	%  RK4
	y_rk4=rk4('de_f',t,y0,h,M+1);
	%e_rk4(im)=norm(y_exact-y_rk4,inf);
	e_rk4(im)=abs(y_exact(M+1)-y_rk4(M+1));

end;

points

% get(gcf)
% set(gcf,'Position', [1203 732 515 307]);

% plot results
loglog(points,e_euler,'--ro','MarkerSize',8)
hold on
loglog(points,e_beuler,'--rd','Color',[0.5 0 0.5],'MarkerSize',8)
loglog(points,e_trap,'--m*','MarkerSize',8)
loglog(points,e_rk4,'--bs','MarkerSize',8)
legend(' Euler',' B Euler',' Trap',' RK4',3);

% define axes used in plot
xlabel('Number of Time Points','FontSize',14,'FontWeight','bold')
ylabel('Error','FontSize',14,'FontWeight','bold')
title(['Logistic Equation'],'FontSize',14,'FontWeight','bold')

% have MATLAB use certain plot options (all are optional)
set(gca,'ytick',[1e-15 1e-12 1e-9 1e-6 1e-3]);
axis([9.99 1e6 1e-15 3e-3])
grid on
set(gca,'MinorGridLineStyle','none')
% Set the fontsize to 14 for the plot
set(gca,'FontSize',14); 
% Set legend font to 14/bold                            		
set(findobj(gcf,'tag','legend'),'FontSize',14,'FontWeight','bold');  

hold off



% right-hand side of DE
function z=de_f(t,y)
r=10;
z=r*y*(1-y);

% exact solution
function y=exact(t,y0,n)
a0=(1-y0)/y0;
r=10;
y=y0;
for i=2:n
	yy=1/(1+a0*exp(-r*t(i)));
	y=[y, yy];
end;

% RK4 method
function ypoints=rk4(f,t,y0,h,n)
y=y0;
ypoints=y0;
for i=2:n
	k1=h*feval(f,t(i-1),y);
	k2=h*feval(f,t(i-1)+0.5*h,y+0.5*k1);
	k3=h*feval(f,t(i-1)+0.5*h,y+0.5*k2);
	k4=h*feval(f,t(i),y+k3);
	yy=y+(k1+2*k2+2*k3+k4)/6;
	ypoints=[ypoints, yy];
	y=yy;
end;

% trap method
function ypoints=trap(f,t,y0,h,n)
%tol=1.0e-7;
tol=10^(-4-n);
y=y0; fold=feval(f,t(1),y);
ypoints=y0;
for i=2:n
	%  secant method
	c=y+0.5*h*fold;
	yb=y;         fb=yb-0.5*h*feval(f,t(i),yb)-c;
	yc=y+0.1*h*fold;  fc=yc-0.5*h*feval(f,t(i),yc)-c;
	while abs(yc-yb) > tol
	   ya=yb; fa=fb;
	   yb=yc; fb=fc;
	   yc=yb-fb*(yb-ya)/(fb-fa);
	   fc=yc-0.5*h*feval(f,t(i),yc)-c;
	end;
	y=yc; fold=feval(f,t(i),y);
	ypoints=[ypoints, y];
end;

% euler method
function ypoints=euler(f,t,y0,h,n)
y=y0;
ypoints=y0;
for i=2:n
	yy=y+h*feval(f,t(i-1),y);
	ypoints=[ypoints, yy];
	y=yy;
end;

% backward euler method
function ypoints=beuler(f,t,y0,h,n)
tol=1.0e-7;
y=y0; fold=feval(f,t(1),y);
ypoints=y0;
for i=2:n
	%  secant method
	yb=y;         fb=yb-h*feval(f,t(i),yb)-y;
	yc=y+2*h*(fold+feval(f,t(i),y+h*fold));  fc=yc-h*feval(f,t(i),yc)-y;
	while abs(yc-yb) > tol
	   ya=yb; fa=fb;
	   yb=yc; fb=fc;
	   yc=yb-fb*(yb-ya)/(fb-fa);
	   fc=yc-h*feval(f,t(i),yc)-y;
	end;
	y=yc; fold=feval(f,t(i),y);
	ypoints=[ypoints, y];
end;